Infiltration.f90 Source File
Implement models to compute infiltration rate
Source Code
!! Implement models to compute infiltration rate !|author: <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a> ! license: <a href="http://www.gnu.org/licenses/">GPL</a> ! !### History ! ! current version 1.1 - 16th June 2016 ! ! | version | date | comment | ! |----------|-------------|----------| ! | 1.0 | 28/May/2014 | Original code | ! | 1.1 | 16/Jun/2016 | added CurveNumber function | ! !### License ! license: GNU GPL <http://www.gnu.org/licenses/> ! !### Module Description ! Implement models to compute infiltration rate: ! ! * [[Philip(function)]] model ! ! * [[GreenAmpt(function)]] model ! ! * [[SCScurveNumber(function)]] model ! MODULE Infiltration USE DataTypeSizes, ONLY : & ! Imported Type Definitions: short, float, double USE LogLib,ONLY : & ! imported routines: Catch USE IniLib, ONLY : & !Imported types: IniList, & !Imported routines: IniOpen, IniClose , & IniReadInt, IniReadReal , & KeyIsPresent, IniReadString, & IniReadDouble, SectionIsPresent USE GridLib, ONLY : & !imported definitions: grid_real, grid_integer, & !imported routines: NewGrid USE StringManipulation, ONLY: & !Imported routines: ToString USE GridOperations, ONLY : & !imported routines: GridByIni, CRSisEqual !, CellArea, & !Imported assignment: !ASSIGNMENT (=) USE SoilProperties, ONLY : & !Imported definitions: SoilType USE DomainProperties, ONLY : & !imported variables: mask USE Richards, ONLY : & !imported routines: SetRichardsBC, SolveRichardsBC, & !imported variables: nsteps, parameters, cell, & dx, S, h0grid, drncum, & infilcum, evapcum, runoffcum IMPLICIT NONE ! Global (i.e. public) Declarations: INTEGER (KIND = short) :: infiltrationModel !infiltration model labels INTEGER (KIND = short), PARAMETER :: SCS_CN = 1 !!Soil Conservation Service Curve Number INTEGER (KIND = short), PARAMETER :: PHILIPEQ = 2 !!Philip's equation INTEGER (KIND = short), PARAMETER :: GREEN_AMPT = 3 !!Green Ampt INTEGER (KIND = short), PARAMETER :: ROSS_BC = 4 !!Richards equations solution by Ross with Broooks and Corey SWRC INTEGER (KIND = short), PARAMETER :: ROSS_VG = 5 !!Richards equations solution by Ross with Van Genuchten SWRC !soil TYPE (SoilType), PRIVATE, ALLOCATABLE :: soils (:) INTEGER (KIND = short) :: soilDivisions !! number of subdivisions of soil layer !common properties TYPE (grid_integer), PRIVATE :: soilTypeMap !read in subroutine SetParametersFromDB TYPE (grid_real), PUBLIC :: ksat !!saturated hydraulic conductivity [m/s] TYPE (grid_real), PUBLIC :: thetar !!residual water content [m3/m3] TYPE (grid_real), PUBLIC :: thetas !!saturated water content [m3/m3] TYPE (grid_real), PUBLIC :: psdi !!pore size distribution index by Brooks and Corey [-] TYPE (grid_real), PUBLIC :: wiltingPoint !!wilting point [m3/m3] TYPE (grid_real), PUBLIC :: fieldCapacity !!field capacity [m3/m3] TYPE (grid_real), PRIVATE :: smax !!maximum soil storage [m] DEBUG non so chi usa questo parametro !used by PHILIPS, ROSS_BC and ROSS_VG TYPE (grid_real), PUBLIC :: psic !!bubbling pressure, air entry value [m] !used by SCS_CN TYPE (grid_real), PUBLIC :: curveNumber !! curve number [-] TYPE (grid_real), PUBLIC :: abstractionRatio !! initial abstraction ratio of the SCS-CN [-] default = 0.2 TYPE (grid_real), PUBLIC :: storativity !! storativity of the SCS-CN [mm], default = 254 TYPE (grid_real), PUBLIC :: sEff !! effective soil retention capacity [mm] TYPE (grid_real), PUBLIC :: raincum !!cumulated gross rainfall [mm] !used by GREEN_AMPT TYPE (grid_real), PUBLIC :: phy !!wetting front soil suction head [m] !used by ROSS_VG TYPE (grid_real), PRIVATE :: ptort !!tortuosity index [-] TYPE (grid_real), PRIVATE :: nvg !!n shape parameter for Van Genuchten SWRC [-] TYPE (grid_real), PRIVATE :: mvg !!m shape parameter for Van Genuchten SWRC [-] TYPE (grid_real), PRIVATE :: ksatMatrix !!saturated conductivity of soil "matrix", !different from ksat when macropores impact. !Implemented only when vG WRC is used !used by PHILIPEQ and GREEN_AMPT TYPE (grid_real), PUBLIC :: cuminf !!cumulative infiltration [m] TYPE (grid_real), PUBLIC :: ism !!soil moisture at the beginning of storm event [m3/m3] !private declarations: INTEGER (KIND = short), PRIVATE :: parameterAssigningMethod !!1 = read each parameter from single file; 2 = assign parameter from soil type map INTEGER (KIND = short), PRIVATE, PARAMETER :: FROM_FILES = 1 INTEGER (KIND = short), PRIVATE, PARAMETER :: FROM_DB = 2 !public routines: PUBLIC :: InfiltrationInit PUBLIC :: Philip PUBLIC :: GreenAmpt PUBLIC :: SCScurveNumber PUBLIC :: Sat2scn !private routines: PRIVATE :: ReadSoilTypes PRIVATE :: SetParametersFromDB PRIVATE :: SetParametersFromFile !======= CONTAINS !======= ! Define procedures contained in this module. !============================================================================== !| Description: ! Initialize evapotranspiration computation SUBROUTINE InfiltrationInit & ! ( inifile, soilMoisture, soilDepth ) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: inifile TYPE (grid_real), INTENT (IN) :: soilMoisture TYPE (grid_real), INTENT (IN) :: soilDepth !local declarations: TYPE (IniList) :: infiltini !--------------------------------------end of declarations--------------------- !open and read configuration file CALL IniOpen (inifile, infiltini) !read infiltration model infiltrationModel = IniReadInt ('model', infiltini) !read number of subdivisions of soil layer soilDivisions = IniReadInt ('ross-divisions', infiltini) !check soil division and infiltration model IF (soilDivisions > 1 .AND. infiltrationModel == SCS_CN .OR. & soilDivisions > 1 .AND. infiltrationModel == PHILIPEQ .OR. & soilDivisions > 1 .AND. infiltrationModel == GREEN_AMPT ) THEN CALL Catch ('warning', 'Infiltration', & 'Infiltration model requires one subdivision of soil layer, divisions forced to 1') soilDivisions = 1 END IF !read parameter assigning method parameterAssigningMethod = IniReadInt ('parameter-assigning-method', infiltini) IF (parameterAssigningMethod == FROM_DB) THEN !check if all necessary info are defined IF (.NOT. KeyIsPresent ('soil-types-file', infiltini) ) THEN CALL Catch ('error', 'Infiltration', & 'missing ''soil-types-file'' keyword in configuration file') END IF IF (.NOT. SectionIsPresent ('soil-type-map', infiltini) ) THEN CALL Catch ('error', 'Infiltration', & 'missing ''soil-type-map'' section in configuration file') END IF !set parameter maps CALL SetParametersFromDB (infiltini, infiltrationModel) ELSE CALL SetParametersFromFile (infiltini, infiltrationModel) END IF !used by Philip model IF (infiltrationModel == PHILIPEQ .OR. & infiltrationModel == GREEN_AMPT ) THEN CALL NewGrid (cuminf, mask, 0.) CALL NewGrid (ism, mask, 0.) END IF !used by SCS-CN IF (infiltrationModel == SCS_CN) THEN CALL NewGrid (raincum, mask, 0.) END IF !when Richards equation solution is required, initialize Ross'parameters IF (infiltrationModel == ROSS_BC) THEN CALL SetRichardsBC (mask, ksat, thetas, thetar, psdi, psic, & soilDepth, soilMoisture, infiltini, soilDivisions) END IF ! Configuration terminated. Deallocate ini database CALL IniClose (infiltini) RETURN END SUBROUTINE InfiltrationInit !============================================================================== !| Description: ! read soil types from external file SUBROUTINE ReadSoilTypes & ! (inifile) IMPLICIT NONE !Arguments with intent(in): CHARACTER (LEN = *), INTENT(IN) :: inifile !!stores configuration information !local declarations TYPE(IniList) :: soilDB INTEGER (KIND = short) :: nsoils !number of soil types to read INTEGER (KIND = short) :: i !------------end of declaration------------------------------------------------ !open and read soil types file CALL IniOpen (inifile, soilDB) !get number of soil types to read IF (KeyIsPresent ('soil-types', soilDB) ) THEN !mandatory key nsoils = IniReadInt ('soil-types', soilDB) ELSE CALL Catch ('error', 'SoilBalance', & 'missing ''soil-types'' section in soil types file') END IF !allocate memory and read data ALLOCATE (soils(nsoils)) DO i = 1, nsoils soils (i) % ksat = IniReadDouble ('ksat', soilDB, section = ToString(i)) soils (i) % thetas = IniReadDouble ('thetas', soilDB, section = ToString(i)) soils (i) % thetar = IniReadDouble ('thetar', soilDB, section = ToString(i)) soils (i) % wp = IniReadDouble ('wp', soilDB, section = ToString(i)) soils (i) % fc = IniReadDouble ('fc', soilDB, section = ToString(i)) soils (i) % psic = IniReadDouble ('psic', soilDB, section = ToString(i)) soils (i) % psdi = IniReadDouble ('psdi', soilDB, section = ToString(i)) soils (i) % phy = IniReadDouble ('phy', soilDB, section = ToString(i)) soils (i) % smax = IniReadDouble ('smax', soilDB, section = ToString(i)) soils (i) % m = IniReadDouble ('m', soilDB, section = ToString(i)) soils (i) % n = IniReadDouble ('n', soilDB, section = ToString(i)) soils (i) % pp = IniReadDouble ('p', soilDB, section = ToString(i)) soils (i) % kx = IniReadDouble ('ksat-matrix', soilDB, section = ToString(i)) soils (i) % c = IniReadDouble ('c', soilDB, section = ToString(i)) soils (i) % s0 = IniReadDouble ('s0', soilDB, section = ToString(i)) soils (i) % cn = IniReadDouble ('cn', soilDB, section = ToString(i)) END DO !deallocate soilDB CALL IniClose (soilDB) RETURN END SUBROUTINE ReadSoilTypes !============================================================================== !| Description: ! set parameter maps from soil database SUBROUTINE SetParametersFromDB & ! (iniDB, model) IMPLICIT NONE ! Arguments with intent(in): TYPE (IniList), INTENT(IN) :: iniDB INTEGER (KIND = short), INTENT(IN) :: model !infiltration model !Local declaration: CHARACTER (LEN = 1000) :: soilTypeFile INTEGER (KIND = short) ::i, j, id !------------end of declaration------------------------------------------------ CALL Catch ('info', 'SoilBalance: ', 'setting soil parameters from database: ', & argument = IniReadString('soil-types-file', iniDB) ) !load soil types soilTypeFile = IniReadString('soil-types-file', iniDB) CALL ReadSoilTypes (soilTypeFile) ! load soil type map CALL GridByIni (iniDB, soilTypeMap, section = 'soil-type-map') !Set parameter maps used by all models !first allocate memory CALL NewGrid (ksat, soilTypeMap) CALL NewGrid (thetar, soilTypeMap) CALL NewGrid (thetas, soilTypeMap) CALL NewGrid (wiltingPoint, soilTypeMap) CALL NewGrid (fieldCapacity, soilTypeMap) CALL NewGrid (psdi, soilTypeMap) !then assigh parameters DO i = 1, soilTypeMap % idim DO j = 1, soilTypeMap % jdim IF (soilTypeMap % mat (i,j) /= soilTypeMap % nodata) THEN id = soilTypeMap % mat (i,j) ksat % mat (i,j) = soils (id) % ksat thetar % mat (i,j) = soils (id) % thetar thetas % mat (i,j) = soils (id) % thetas wiltingPoint % mat (i,j) = soils (id) % wp fieldCapacity % mat (i,j) = soils (id) % fc psdi % mat (i,j) = soils (id) % psdi END IF END DO END DO IF (model == SCS_CN) THEN !read supplementary parameters required by Curve Number !first allocate memory CALL NewGrid (curveNumber, soilTypeMap) CALL NewGrid (abstractionRatio, soilTypeMap) CALL NewGrid (storativity, soilTypeMap) !then assign parameters DO i = 1, soilTypeMap % idim DO j = 1, soilTypeMap % jdim IF (soilTypeMap % mat (i,j) /= soilTypeMap % nodata) THEN id = soilTypeMap % mat (i,j) curveNumber % mat (i,j) = soils (id) % cn abstractionRatio % mat (i,j) = soils (id) % c storativity % mat (i,j) = soils (id) % s0 END IF END DO END DO END IF IF (model == PHILIPEQ) THEN !read supplementary parameters required by Philips !first allocate memory CALL NewGrid (psic, soilTypeMap) !then assigh parameters DO i = 1, soilTypeMap % idim DO j = 1, soilTypeMap % jdim IF (soilTypeMap % mat (i,j) /= soilTypeMap % nodata) THEN id = soilTypeMap % mat (i,j) psic % mat (i,j) = soils (id) % psic END IF END DO END DO END IF IF (model == GREEN_AMPT) THEN !read supplementary parameters required by Green Ampt !first allocate memory CALL NewGrid (phy, soilTypeMap) !then assigh parameters DO i = 1, soilTypeMap % idim DO j = 1, soilTypeMap % jdim IF (soilTypeMap % mat (i,j) /= soilTypeMap % nodata) THEN id = soilTypeMap % mat (i,j) phy % mat (i,j) = soils (id) % phy END IF END DO END DO END IF IF (model == ROSS_BC) THEN !read supplementary parameters required by Ross Brooks and Corey !first allocate memory CALL NewGrid (psic, soilTypeMap) !then assigh parameters DO i = 1, soilTypeMap % idim DO j = 1, soilTypeMap % jdim IF (soilTypeMap % mat (i,j) /= soilTypeMap % nodata) THEN id = soilTypeMap % mat (i,j) psic % mat (i,j) = soils (id) % psic END IF END DO END DO END IF IF (model == ROSS_VG) THEN !read supplementary parameters required by Ross Van Genuchten !first allocate memory CALL NewGrid (psic, soilTypeMap) CALL NewGrid (nvg, soilTypeMap) CALL NewGrid (mvg, soilTypeMap) CALL NewGrid (ptort, soilTypeMap) CALL NewGrid (ksatMatrix, soilTypeMap) !then assigh parameters DO i = 1, soilTypeMap % idim DO j = 1, soilTypeMap % jdim IF (soilTypeMap % mat (i,j) /= soilTypeMap % nodata) THEN id = soilTypeMap % mat (i,j) psic % mat (i,j) = soils (id) % psic nvg % mat (i,j) = soils (id) % n mvg % mat (i,j) = soils (id) % m mvg % mat (i,j) = soils (id) % m ksatMatrix % mat (i,j) = soils (id) % kx END IF END DO END DO END IF RETURN END SUBROUTINE SetParametersFromDB !============================================================================== !| Description: ! set parameter maps from separate files SUBROUTINE SetParametersFromFile & ! (iniDB, model) IMPLICIT NONE ! Arguments with intent(in): TYPE (IniList), INTENT(IN) :: iniDB INTEGER (KIND = short), INTENT(IN) :: model !infiltration model !Local declaration: !------------end of declaration------------------------------------------------ !Set parameter maps used by all models !Saturated Conductivity IF (SectionIsPresent('conductivity', iniDB)) THEN CALL GridByIni (iniDB, ksat, section = 'conductivity') CALL Catch ('info', 'SoilBalance: ', 'loading saturated conductivity from file: ', & argument = IniReadString('file', iniDB, section = 'conductivity')) IF ( .NOT. CRSisEqual (mask = mask, grid = ksat, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in conductivity map' ) END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing conductivity section in configuration file' ) END IF !Saturated water content IF (SectionIsPresent('saturated-water-content', iniDB)) THEN CALL GridByIni (iniDB, thetas, section = 'saturated-water-content') CALL Catch ('info', 'SoilBalance: ', 'loading saturated water content from file: ', & argument = IniReadString('file', iniDB, section = 'saturated-water-content')) IF ( .NOT. CRSisEqual (mask = mask, grid = thetas, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in saturated-water-content' ) END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing saturated-water-content section in configuration file' ) END IF !Residual water content IF (SectionIsPresent('residual-water-content', iniDB)) THEN CALL GridByIni (iniDB, thetar, section = 'residual-water-content') CALL Catch ('info', 'SoilrBalance: ', 'loading residual water content from file: ', & argument = IniReadString('file', iniDB, section = 'residual-water-content')) IF ( .NOT. CRSisEqual (mask = mask, grid = thetar, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in residual-water-content' ) END IF ELSE CALL Catch ('error', 'SoilWaterBalance: ', & 'missing residual-water-content section in configuration file' ) END IF !Wilting Point IF (SectionIsPresent('wilting-point', iniDB)) THEN CALL GridByIni (iniDB, wiltingPoint, section = 'wilting-point') CALL Catch ('info', 'SoilBalance: ', 'loading wilting point from file: ', & argument = IniReadString('file', iniDB, section = 'wilting-point')) IF ( .NOT. CRSisEqual (mask = mask, grid = wiltingPoint, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in wilting-point' ) END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing wilting-point section in configuration file' ) END IF !Field Capacity IF (SectionIsPresent('field-capacity', iniDB)) THEN CALL GridByIni (iniDB, fieldCapacity, section = 'field-capacity') CALL Catch ('info', 'SoilBalance: ', 'loading field capacity from file: ', & argument = IniReadString('file', iniDB, section = 'field-capacity')) IF ( .NOT. CRSisEqual (mask = mask, grid = fieldCapacity, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in field-capacity' ) END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing field-capacity section in configuration file' ) END IF !Pore size distribution index IF (SectionIsPresent('pore-size-index', iniDB)) THEN CALL GridByIni (iniDB, psdi, section = 'pore-size-index') CALL Catch ('info', 'SoilBalance: ', 'loading pore size index from file: ', & argument = IniReadString('file', iniDB, section = 'pore-size-index')) IF ( .NOT. CRSisEqual (mask = mask, grid = psdi, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in pore-size-index') END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing pore-size-idex section in configuration file' ) END IF !bubbling pressure, air entry value IF (SectionIsPresent('bubble-pressure', iniDB)) THEN CALL GridByIni (iniDB, psic, section = 'bubble-pressure') CALL Catch ('info', 'SoilBalance: ', 'loading bubbling pressure from file: ', & argument = IniReadString('file', iniDB, section = 'bubble-pressure')) IF ( .NOT. CRSisEqual (mask = mask, grid = psic, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in bubbling pressure') END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing bubble-pressure section in configuration file' ) END IF IF (model == SCS_CN) THEN !read supplementary parameters required by Curve Number !Curve Number IF (SectionIsPresent('curve-number', iniDB)) THEN CALL GridByIni (iniDB, curveNumber, section = 'curve-number') CALL Catch ('info', 'SoilBalance: ', 'loading curve number from file: ', & argument = IniReadString('file', iniDB, section = 'curve-number')) IF ( .NOT. CRSisEqual (mask = mask, grid = curveNumber, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in curve number') END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing curve-number section in configuration file' ) END IF !Initial abstraction ratio IF (SectionIsPresent('abstraction-ratio', iniDB)) THEN CALL GridByIni (iniDB, abstractionRatio, section = 'abstraction-ratio') CALL Catch ('info', 'SoilBalance: ', 'loading initial abstraction ratio from file: ', & argument = IniReadString('file', iniDB, section = 'abstraction-ratio')) IF ( .NOT. CRSisEqual (mask = mask, grid = abstractionRatio, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in initial abstraction ratio') END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing abstraction-ratio section in configuration file' ) END IF !SCS-CN storativity IF (SectionIsPresent('storativity', iniDB)) THEN CALL GridByIni (iniDB, storativity, section = 'storativity') CALL Catch ('info', 'SoilBalance: ', 'loading storativity from file: ', & argument = IniReadString('file', iniDB, section = 'storativity')) IF ( .NOT. CRSisEqual (mask = mask, grid = storativity, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in storativity') END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing storativity section in configuration file' ) END IF END IF !IF (model == PHILIPEQ) THEN !read supplementary parameters required by Philips ! ! !bubbling pressure, air entry value ! IF (SectionIsPresent('bubble-pressure', iniDB)) THEN ! CALL GridByIni (iniDB, psic, section = 'bubble-pressure') ! CALL Catch ('info', 'SoilBalance: ', 'loading bubbling pressure from file: ', & ! argument = IniReadString('file', iniDB, section = 'bubble-pressure')) ! IF ( .NOT. CRSisEqual (mask = mask, grid = psic, checkCells = .TRUE.) ) THEN ! CALL Catch ('error', 'SoilBalance', & ! 'wrong spatial reference in bubbling pressure') ! END IF ! ELSE ! CALL Catch ('error', 'SoilBalance: ', & ! 'missing bubble-pressure section in configuration file' ) ! END IF ! !END IF IF (model == GREEN_AMPT) THEN !read supplementary parameters required by Green Ampt !wetting front soil suction head IF (SectionIsPresent('front-suction-head', iniDB)) THEN CALL GridByIni (iniDB, phy, section = 'front-suction-head') CALL Catch ('info', 'SoilBalance: ', 'loading wetting front suction head from file: ', & argument = IniReadString('file', iniDB, section = 'front-suction-head')) IF ( .NOT. CRSisEqual (mask = mask, grid = phy, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in wetting front suction head') END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing front-suction-head section in configuration file' ) END IF END IF IF (model == ROSS_BC) THEN !read supplementary parameters required by Ross Brooks and Corey !bubbling pressure, air entry value IF (SectionIsPresent('bubble-pressure', iniDB)) THEN CALL GridByIni (iniDB, psic, section = 'bubble-pressure') CALL Catch ('info', 'SoilBalance: ', 'loading bubbling pressure from file: ', & argument = IniReadString('file', iniDB, section = 'bubble-pressure')) IF ( .NOT. CRSisEqual (mask = mask, grid = psic, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in bubbling pressure') END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing bubble-pressure section in configuration file' ) END IF END IF IF (model == ROSS_VG) THEN !read supplementary parameters required by Ross Van Genuchten !bubbling pressure, air entry value IF (SectionIsPresent('bubble-pressure', iniDB)) THEN CALL GridByIni (iniDB, psic, section = 'bubble-pressure') CALL Catch ('info', 'SoilBalance: ', 'loading bubbling pressure from file: ', & argument = IniReadString('file', iniDB, section = 'bubble-pressure')) IF ( .NOT. CRSisEqual (mask = mask, grid = psic, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in bubbling pressure') END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing bubble-pressure section in configuration file' ) END IF !n shape parameter for Van Genuchten SWRC IF (SectionIsPresent('n-van-genuchten', iniDB)) THEN CALL GridByIni (iniDB, nvg, section = 'n-van-genuchten') CALL Catch ('info', 'SoilBalance: ', 'loading n shape parameter from file: ', & argument = IniReadString('file', iniDB, section = 'n-van-genuchten')) IF ( .NOT. CRSisEqual (mask = mask, grid = nvg, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in n shape parameter') END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing n-van-genuchten section in configuration file' ) END IF !m shape parameter for Van Genuchten SWRC IF (SectionIsPresent('m-van-genuchten', iniDB)) THEN CALL GridByIni (iniDB, mvg, section = 'm-van-genuchten') CALL Catch ('info', 'SoilBalance: ', 'loading m shape parameter from file: ', & argument = IniReadString('file', iniDB, section = 'm-van-genuchten')) IF ( .NOT. CRSisEqual (mask = mask, grid = mvg, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in m shape parameter') END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing m-van-genuchten section in configuration file' ) END IF !tortuosity index IF (SectionIsPresent('tortuosity-index', iniDB)) THEN CALL GridByIni (iniDB, ptort, section = 'tortuosity-index') CALL Catch ('info', 'SoilBalance: ', 'loading tortuosity index from file: ', & argument = IniReadString('file', iniDB, section = 'tortuosity-index')) IF ( .NOT. CRSisEqual (mask = mask, grid = ptort, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in tortuosity index') END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing tortuosity-index section in configuration file' ) END IF !saturated conductivity of soil "matrix" IF (SectionIsPresent('conductivity-matrix', iniDB)) THEN CALL GridByIni (iniDB, ksatMatrix, section = 'conductivity-matrix') CALL Catch ('info', 'SoilBalance: ', 'loading soil matrix conductivity from file: ', & argument = IniReadString('file', iniDB, section = 'conductivity-matrix')) IF ( .NOT. CRSisEqual (mask = mask, grid = ksatMatrix, checkCells = .TRUE.) ) THEN CALL Catch ('error', 'SoilBalance', & 'wrong spatial reference in soil matrix conductivity') END IF ELSE CALL Catch ('error', 'SoilBalance: ', & 'missing conductivity-matrix section in configuration file' ) END IF END IF RETURN END SUBROUTINE SetParametersFromFile !============================================================================== !| calculates the actual rate of infiltration (m/s) according to ! Curve Number model modified for continuous simulation ! ! References: ! ! Ravazzani, G., Mancini, M., Giudici, I., & Amadio, P.. (2007). ! Effects of soil moisture parameterization on a real- time flood ! forecasting system based on rainfall thresholds. ! IAHS publication, 313, 407–416. ! ! Rabuffetti, D., Ravazzani, G., Corbari, C., & Mancini, M.. (2008). ! Verification of operational quantitative discharge forecast (QDF) ! for a regional warning system – the AMPHORE case studies in the ! upper Po river. Natural hazards and earth system sciences, 8, 161–173. ! ! Ravazzani, G., Amengual, A., Ceppi, A., Homar, V., Romero, R., ! Lombardi, G., & Mancini, M.. (2016). Potentialities of ensemble ! strategies for flood forecasting over the Milano urban area. ! Journal of hydrology, 539, 237-253. ! FUNCTION SCScurveNumber & ! (rain, raincum, c, s, dt, runoff) & ! RESULT (inf) IMPLICIT NONE !Arguments with intent in REAL (KIND = float), INTENT(IN) :: rain !!rainfall rate (m/s) REAL (KIND = float), INTENT(IN) :: c !!initial abstraction ratio (-) REAL (KIND = float), INTENT(IN) :: s !!soil retention capacity (mm) INTEGER (KIND = short), INTENT(IN) :: dt !!time step (s) !Arguments with intent inout: REAL (KIND = float), INTENT(INOUT) :: raincum !!cumulated rainfall from start of storm (mm) !! is used and updated to be ready for the !! following step !Arguments with intent out REAL (KIND = double), INTENT(OUT) :: runoff !!runoff rate (m/s) !local declarations: REAL (KIND = float) :: inf !!infiltration rate (m/s) REAL (KIND = float) :: raincump !!cumulative precipitation of previous time step (mm) REAL (KIND = float) :: runoffp !!cumulative runoff of previous time step (mm) !-------------------------end of declarations---------------------------------- !if retention capacity = 0 (soil is saturated) ! rainfall is transformed to runoff IF (s == 0.) THEN raincum = raincum + rain * 1000. * dt !update cumulative rainfall runoff = rain inf = 0. RETURN !no need to continue END IF !calculate runoff and infiltration rate !update cumulative precipitation raincump = raincum !save cumulated precipitation amount of previous time step raincum = raincum + rain * dt * 1000. !(mm) !runoff at current time IF(raincum >= c*s) THEN runoff = ((raincum - c*s )**2.) / (raincum + (1.-c) * s) ELSE runoff = 0. END IF !runoff at previous time IF(raincump >= c*S) THEN runoffp = ((raincump - c*S )**2.) / (raincump + (1.-c) * s) ELSE runoffp = 0. END IF !actual runoff in m/s runoff = (runoff - runoffp) / 1000. / dt !actual infiltration rate in m/s inf = rain - runoff RETURN END FUNCTION SCScurveNumber !============================================================================== !| compute actual soil retention capacity S of the SCS-CN method ! given soil saturation SUBROUTINE Sat2scn & ! (cn2, s0, soilsat, scn) IMPLICIT NONE !arguments with intent in: REAL (KIND = float), INTENT(IN) :: cn2 !!curve number for mean soil moisture REAL (KIND = float), INTENT(IN) :: s0 !!storativity (mm) REAL (KIND = float), INTENT(IN) :: soilsat !!soil saturation !arguments with intent out: REAL (KIND = float), INTENT(OUT) :: scn !!soil retention capacity [mm] !local declarations: REAL (KIND = float) :: cn1, cn3, s1, s2, s3 !-------------------------end of declarations---------------------------------- !compute cn1 and cn3 from cn2 cn1 = cn2 - (20. * (100. - cn2) ) / & (100. - cn2 + EXP(2.533 - .0636 * (100. - cn2) ) ) IF (cn1 < 0.4 * CN2) cn1 = .4 * cn2 cn3 = cn2 * EXP( .00673 * (100. - cn2) ) !compute s1, s2, and s3 s1 = s0 * (100. / cn1 - 1.) s2 = s0 * (100. / cn2 - 1.) s3 = s0 * (100. / cn3 - 1.) !compute actual S scn = s1 - soilsat * (s1 - s3) RETURN END SUBROUTINE sat2scn !============================================================================== !| calculates the actual rate of infiltration (m/s) according to ! Philip's equation. Use time compression approximation. ! ! References: ! ! Philip, J. R. 1954. An infiltration equation with physical ! significance. Soil Science, 77: 153-157. ! ! Milly, P. C. D., An event-based simulation model of moisture and ! energy fluxes at a bare soil surface, ! Water Resour. Res., 22, 1680-! 692, 1986 ! ! Sivapalan, M., Beven, K., and Wood, E. F. (1987), On hydrologic similarity: ! 2 . A scaled model of storm runoff production, Water Resour. ! Res., 23( 12), 2266– 2278, doi:10.1029/WR023i012p02266. ! FUNCTION Philip & ! (rain, psic, ksat, theta, thetas, thetar, psdi, itheta, dt, cuminf) & ! RESULT (inf) IMPLICIT NONE !Arguments with intent in REAL (KIND = float), INTENT(IN) :: rain !!rainfall rate (m/s) REAL (KIND = float), INTENT(IN) :: psic !!bubbling pressure (m) REAL (KIND = float), INTENT(IN) :: ksat !!saturated hydraulic conductivity (m/s) REAL (KIND = float), INTENT(IN) :: theta !!volumetric water content (m3/m3) REAL (KIND = float), INTENT(IN) :: thetas !!saturated volumetric water content (m3/m3) REAL (KIND = float), INTENT(IN) :: thetar !!residual volumetric water content (m3/m3) REAL (KIND = float), INTENT(IN) :: psdi !!Brooks & Corey pore size distribution index (-) INTEGER (KIND = short), INTENT(IN) :: dt !!time step (s) REAL (KIND = float), INTENT(IN) :: itheta !!initial soil moisture !Arguments with intent inout: REAL (KIND = float), INTENT(INOUT) :: cuminf !!cumulative infiltration (m) !local declarations: REAL (KIND = float) :: inf !!infiltration rate (m/s) REAL (KIND = float) :: sorp !!sorptivity REAL (KIND = float) :: cc = 0. !!gravity term in Philip's infiltration equation REAL (KIND = float) :: deltrz !! difference between soil moisture at beginning !! of storm event and residual soil moisture REAL (KIND = float) :: bcgamm !!2.0 + 3.0 * psdi REAL (KIND = float) :: infcap !!infiltration capacity (m/s) !-------------------------end of declarations---------------------------------- ! Sivapalan et al. (1987) bcgamm = 2.0 + 3.0 * psdi sorp = (((2. * ksat *( (thetas-itheta)**2.) * psic ) / & (thetas-thetar))*((1./(bcgamm + 0.5 * psdi - 1.)) + & ((thetas-thetar)/ (thetas-itheta) ) ) ) ** 0.5 deltrz = itheta - thetar IF (deltrz <= 0.) deltrz = 0. cc = 0.5 * (1. + ( (deltrz / (thetas - thetar))**(bcgamm / psdi) ) ) !If soil is saturated then set infiltration to zero IF (theta >= thetas) THEN inf = 0. !If surface is unsaturated then calculate infiltration ELSE !If first step with infiltration then all precipitation is !infiltrated to avoid division by zero IF (cuminf == 0.) THEN inf = rain !calculate infiltration capacity as a function of cumulative ! infiltration for all other time steps of storm (Milly, 1986) ELSE infcap = cc * ksat * (1. + (1./ ((( 1. + (( 4. * cc * ksat * cuminf ) / & (sorp**2.)))**0.5) - 1.))) !Take the actual infiltration rate as the minimum of the ! precipitation rate or the infiltration capacity. inf = MIN (infcap,rain) END IF !update cumulative infiltration cuminf = cuminf + inf * dt END IF RETURN END FUNCTION Philip !============================================================================== !| calculates the actual rate of infiltration (m/s) according to ! Green-Ampt equation for unsteady rainfall. ! ! References: ! ! Green, WH., and Ampt, GA., 1911. Studies on soil physics I. ! Flow of air and water through soils. J. Agric. Sci. 4: 1- 24. ! ! Chow, V., Maidment, D., Mays, L. (1988) Applied Hydrology, ! Surface water, pp. 142-145, McGraw-Hill ! FUNCTION GreenAmpt & ! (storm, rain, ksat, theta, thetas, thetar, phy, itheta, dt, cuminf) & ! RESULT (inf) IMPLICIT NONE !Arguments with intent in INTEGER (KIND = short), INTENT(IN) :: storm !!when = 1 first steo of storm event REAL (KIND = float), INTENT(IN) :: rain !!rainfall rate (m/s) REAL (KIND = float), INTENT(IN) :: ksat !!saturated hydraulic conductivity (m/s) REAL (KIND = float), INTENT(IN) :: theta !!volumetric water content (m3/m3) REAL (KIND = float), INTENT(IN) :: thetas !!saturated volumetric water content (m3/m3) REAL (KIND = float), INTENT(IN) :: thetar !!residual volumetric water content (m3/m3) REAL (KIND = float), INTENT(IN) :: phy !!suction head across the wetting front (m) INTEGER (KIND = short), INTENT(IN) :: dt !!time step (s) REAL (KIND = float), INTENT(IN) :: itheta !!initial soil moisture !Arguments with intent inout: REAL (KIND = float), INTENT(INOUT) :: cuminf !!cumulative infiltration (m) !local declarations: REAL (KIND = float) :: inf !!infiltration rate (m/s) REAL (KIND = float) :: deltrz !! difference between soil moisture at beginning !! of storm event and residual soil moisture REAL (KIND = float) :: test !!test value for the cumulative infiltration (m) REAL (KIND = float) :: ftest !!test value for the cumulative infiltration (m) INTEGER (KIND = short) :: i !loop INTEGER (KIND = short) :: maxiter = 100 !maximum number of iteration for ! computing cumulative infiltration REAL (KIND = float) :: tol = 0.00001 !tolerance to compute cumulative infiltration (m) REAL (KIND = float) :: deltatfirst !!for ponding time (case 3) REAL (KIND = float) :: Fponding !!cumulative infiltration at ponding (case 3) REAL (KIND = float) :: deltatpond !! new deltat calculated for ponding (case 3) !-------------------------end of declarations---------------------------------- !if this is the first step of storm event: ! initialize variables IF ( storm == 1) THEN deltrz = thetas - theta IF (deltrz <= 0.) THEN deltrz = 0. END IF !assume all precipitation is infiltrated cuminf = rain * dt inf = rain RETURN END IF !If soil is saturated then set infiltration to zero IF (theta >= thetas) THEN inf = 0. RETURN !If surface is unsaturated then calculate infiltration ELSE IF (storm == 2) THEN !storm is running !calculate potential infiltration rate at the beginning of the interval !from the known value of cumulative infiltration (eq 5.4.1) inf = ksat * (phy * (thetas-itheta) / cuminf + 1.) !compare potential infiltration to rainfall rate IF (inf <= rain) THEN !ponding occurs throughout the interval (case 1) test = Ksat*dt !first guess for cumulative infiltration DO i = 1, maxiter IF (i == maxiter) THEN CALL Catch ('warning', 'Infiltration', 'maximum number of iterations & reached while updating cumulative infiltration & with green-ampt case 1' ) END IF ftest = cuminf + ksat*dt + (phy*(thetas-itheta))* & LOG ((test+(phy*(thetas-itheta)))/(cuminf+(phy*(thetas-itheta)))) IF (ABS(ftest-test) <= tol) THEN cuminf = ftest !compute infiltration rate from cumulative infiltration inf = ksat * (phy * (thetas-itheta) / cuminf + 1.) IF (inf > rain) inf = rain EXIT ELSE test = ftest END IF END DO ELSE !no ponding at the beginning of the interval !calculate tentative values Assuming that all the rain !water is penetrated into the soil test = cuminf + rain * dt inf = ksat * (phy * (thetas-itheta) / test + 1.) IF (inf <= rain) THEN !ponding occurs during the interval (case 3) Fponding = ksat * phy * (thetas-itheta) / (rain - ksat) deltatfirst = (Fponding - cuminf) / rain deltatpond = dt - deltatfirst !compute cumulative infiltration by substituting cuminf = Fponding !and dt = deltatpond test = Fponding !Ksat*dt !first guess for cumulative infiltration DO i = 1, maxiter IF (i == maxiter) THEN CALL Catch ('warning', 'Infiltration', 'maximum number of iterations & reached while updating cumulative infiltration & with green-ampt case 3' ) END IF ftest = Fponding + ksat*deltatpond + (phy*(thetas-itheta))* & LOG ((test+(phy*(thetas-itheta)))/(Fponding+(phy*(thetas-itheta)))) IF (ABS(ftest-test) <= tol) THEN cuminf = ftest !compute infiltration rate from cumulative infiltration inf = ksat * (phy * (thetas-itheta) / cuminf + 1.) IF (inf > rain) inf = rain EXIT ELSE test = ftest END IF END DO ELSE !no ponding throghout the interval (case 2) cuminf = test inf = rain END IF END IF ELSE inf = 0. END IF END IF END FUNCTION GreenAmpt END MODULE Infiltration